There are four species of Anax (Emperor dragonflies) found in South Africa. Their ranges are extensive. This is a handy bit of code to visualize their distributions using R and GIS wizardry.
First, let’s set up our work space. These are the packages that we’ll need to run this code, check that they’re installed and in your library:
my_packages <- c("rmarkdown", "tidyverse", "sf", "rinat", "raster", "lwgeom", "leaflet","leafpop","rosm", "ggspatial", "prettymapr", "mapview", "rnaturalearth","readxl","htmltools","bslib")
lapply(my_packages, require, character.only = TRUE)
We’ll need some data. iNaturalist is a great place to find
observational data on organisms around the world, including dragonflies.
The data for this exercise was downloaded and stored in the GitHub
repository - you can find it here.
The four species of dragonflies are:
| Common Name | Scientific Name | |
|---|---|---|
| Blue Emperor | Anax imperator | |
| Orange Emperor | Anax speratus | |
| Vagrant Emperor | Anax ephippiger | |
| Black Emperor | Anax tristis |
The species names will be used for the respective data sets
below:
We’ll begin by reading in the data:
imperator <- read.csv("iNatData/obs_imperator.csv")
speratus <- read.csv("iNatData/obs_speratus.csv")
ephippiger <-read.csv("iNatData/obs_ephippiger.csv")
tristis <- read.csv("iNatData/obs_tristis.csv")
Then we’ll clean it up. Because iNaturalist makes use of citizen
scientists to log their observations, the data isn’t always accurate or
complete.
We’re filtering the data frames to remove observations without
locations, and those with positional accuracy grater than 500m
(dragonflies move around freely, so I’ve allowed for quite a bit
inaccuracy here). We’ll also only include ‘research grade’ data, which
has been somewhat verified by other users.
imperator <- imperator %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
speratus <- speratus %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
ephippiger <- ephippiger %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
tristis <- tristis %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
In order for the GIS wizardry to work, we’ll convert these
tabular data frames to spatial objects using the Sf package:
imperatorSf <- st_as_sf(imperator, coords = c("longitude", "latitude"), crs = 4326)
speratusSf <- st_as_sf(speratus, coords = c("longitude", "latitude"), crs = 4326)
ephippigerSf <- st_as_sf(ephippiger, coords = c("longitude", "latitude"), crs = 4326)
tristisSf <- st_as_sf(tristis, coords = c("longitude", "latitude"), crs = 4326)
Then we’ll use the following code to check that the data is
indeed of the Sf class, and that the correct co-ordinate reference
system is used for each:
cat("<details><summary>Click to expand the meaty results of these checks</summary>\n")
class(imperatorSf)
[1] “sf” “data.frame”
class(speratusSf)
[1] “sf” “data.frame”
class(ephippigerSf)
[1] “sf” “data.frame”
class(tristisSf)
[1] “sf” “data.frame”
st_crs(imperatorSf)
Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]
st_crs(speratusSf)
Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]
st_crs(ephippigerSf)
Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]
st_crs(tristisSf)
Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]
cat("</details>\n")
Now that we’re confident that our data are accurate, research
grade, and spatial objects, we can visualize them on a map:
ggplot() +
annotation_map_tile(type = "osm", progress = "none") +
geom_sf(data = imperatorSf, color = "#83aed1") +
ggtitle(expression("Distribution of " * italic("Anax") * italic(" imperator")))
ggplot() +
annotation_map_tile(type = "osm", progress = "none") +
geom_sf(data = speratusSf, color = "#e1a56d") +
ggtitle(expression("Distribution of " * italic("Anax") * italic("speratus")))
ggplot() +
annotation_map_tile(type = "osm", progress = "none") +
geom_sf(data = ephippigerSf, color = "#afbd6c") +
ggtitle(expression("Distribution of " * italic("Anax") * italic(" ephippiger")))
ggplot() +
annotation_map_tile(type = "osm", progress = "none") +
geom_sf(data = tristisSf, color = "#9980ab") +
ggtitle(expression("Distribution of " * italic("Anax") * italic(" tristis")))
Great! These maps show the distribution of these species well, and
the colour coding helps to differentiate them. However, these maps
aren’t actually very useful to look at the observational data of one
species in relation to another, nor does it allow us to visulise the
data on any meaningful scale.
An interactive map, with all the
species on it, would be much more useful! So we’ll use the Leaflet
package to create a better map:
leaflet() %>%
addTiles(group = "Default") %>%
addCircleMarkers(data = imperatorSf,
group = "Anax imperator",
radius = 2,
color = "#83aed1") %>%
addCircleMarkers(data = speratusSf,
group = "Anax speratus",
radius = 2,
color = "#e1a56d") %>%
addCircleMarkers(data = ephippigerSf,
group = "Anax ephippiger",
radius = 2,
color = "#afbd6c") %>%
addCircleMarkers(data = tristisSf,
group = "Anax tristis",
radius = 2,
color = "#9980ab") %>%
addLegend(position = "topright",
colors = c("#83aed1", "#e1a56d", "#afbd6c", "#9980ab"),
labels = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"))
Much better!
Although, it would be easier to visualize this
data if we could choose which species to display. So lets do that again,
adding a layer control to choose which species’ data we’d like to see:
And while we’re at it, it would make sense to see more details for each
observation - perhaps even a clickable link to the iNaturalist page?
This chunk of code will create clickable links for us to add to
popup labels:
limperatorSf <- imperatorSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>"))
lsperatusSf <- speratusSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>"))
lephippigerSf <- ephippigerSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>"))
ltristisSf <- tristisSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>"))
And now we’re ready to create our user-friendly, and very
pretty, interactive map:
leaflet() %>%
addTiles(group = "Default") %>%
addCircleMarkers(data = imperatorSf,
group = "Blue Emperor",
radius = 2,
color = "#83aed1",
popup = popupTable(limperatorSf,zcol = c("click_url"))) %>%
addCircleMarkers(data = speratusSf,
group = "Orange Emperor",
radius = 2,
color = "#e1a56d",
popup = popupTable(lsperatusSf,zcol = c("click_url"))) %>%
addCircleMarkers(data = ephippigerSf,
group = "Vagrant Emperor",
radius = 2,
color = "#afbd6c",
popup = popupTable(lephippigerSf,zcol = c("click_url"))) %>%
addCircleMarkers(data = tristisSf,
group = "Black Emperor",
radius = 2,
color = "#9980ab",
popup = popupTable(ltristisSf,zcol = c("click_url"))) %>%
addLayersControl( overlayGroups = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"),
options = layersControlOptions(collapsed = FALSE) ) %>%
addLegend(position = "topright",
colors = c("#83aed1", "#e1a56d", "#afbd6c", "#9980ab"),
labels = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"))
It’s amazing that such small creatures have such a large range!
There are at least 9 species of Anax dragonflies that are
scientifically documented to migrate. While we still have much to learn
about these impressive journeys, this distribution map shows that it’s
likely that at least some populations of these four species migrate too.
A special mention to everyone who uses iNaturalist, citizen science
makes observational data projects like this possible! And thanks to Jasper Slingsby for showing me how to
make this, and for all his patience.